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SUMMARY 

We present a multigrid solver for the exponential fitting method, applied to the current con- 
tinuity equations of semiconductor device simulation in two dimensions. The exponential fitting 
method is based on a mixed finite element discretization using the lowest-order Raviart-Thomas 
triangular element. This discretization method yields a good approximation of front layers and 
guarantees current conservation. The corresponding stiffness matrix is an M- matrix. ” Standard’ 
multigrid solvers, however, cannot be applied to the resulting system, as this is dominated by 
an unsymmetric part, which is due to the presence of strong convection in part of the domain. 

To overcome this difficulty, we explore the connection between Raviart-Thomas mixed methods 
and the nonconforming Crouzeix-Raviart finite element discretization. In this way we can con- 
struct nonstandard prolongation and restriction operators using easily computable weighted L 2 - 
projections based on suitable quadrature rules and the upwind effects of the discretization. The 
resulting multigrid algorithm shows very good results, even for real-world problems and for lo- 
cally refined grids. 

1. INTRODUCTION 

The exponential fitting method applied to the current continuity equations is based on a 
mixed finite element discretization using the lowest-order Raviart-Thomas triangular element 
[1], This discretization yields a good approximation of front layers and guarantees current con- 
servation. The corresponding scheme results in a large sparse system of equations, which is dom- 
inated by an unsymmetric part. When applying multigrid algorithms to the resulting system (7), 
the most difficult part is the construction of suitable prolongation and restriction operators. Us- 
ing the connection between Raviart-Thomas mixed methods and the nonconforming Crouzeix- 
Raviart finite element discretization, we overcome this difficulty. 

In § 2 we give some results from [2] concerning the mixed finite element discretization. We 
determine the resulting system and show the interrelation with a nonconforming finite element 
method. §3 deals with the solution of the system of linear equations by our multigrid solver. 
First we construct easily computable L 2 — projections, based on suitable quadrature rules and 
the upwind effect of the discretization. Due to the presence of strong convection in part of the 
domain it is also necessary to consider special smoothers for the multigrid algorithm. We use a 
minimal residual method with ILU preconditioning. The results of the numerical tests are given 
in §4. 
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2. THE EXPONENTIAL FITTING METHOD FOR CURRENT CONTINUITY EQUATIONS 


2.1. Mixed Finite Element Method 


Let Q C 1R 2 be a connected, bounded and polygonal domain. for m G N, and 

L 2 (Q) := H°(Q) denote the usual Sobolev and Lebesque spaces equipped with the norm 



For / G L 2 (Q) and g € L 2 (r 0 ), To C 50 closed with positive length, we consider the current 
continuity equation, as given in [3]: 

Find u G H l (Sl) such that 


div (grad u + u grad ip) = f 

u = g 

du dip _ 


in O C IR 2 , 
on r 0 C 50, 

on Ti = 50 \ Tq. 


( 1 ) 


The current is defined by 3 - grad u + u grad ip. Here, ip G H 1 ( O) is a given bounded function. 
To discretize problem (1) we introduce the classical method of changing variables from u to the 
socalled Slotboom variable p [3] 

p = u. 


T his results in the following symmetric form of problem (1): 


Find p G H 1 { O) such that 


div (e ^grad p) = f 


X := 


e^g 


2 - 


in O C IR 2 , 
on r 0 C 50, 

on Tj = 50 \ r 0 . 


( 2 ) 


Let {7*;}fc>o be a regular sequence of decompositions of O into triangles. Denote by hk the 
longest side of all triangles T G 7*,. The set of edges of 7*. is denoted by £*, where are the 
boundary edges and ££ = £ k \ £% are all interelement boundaries. Denote by m e the midpoint of 
an edge e of £*. Moreover, let P m , m > 0, be the space of all polynomials of degree not greater 
than m. Following [1], we use the lowest order Raviart-Thomas mixed finite element to discretize 
(2). Therefore we define the following set of polynomial vectors 

RT(T) := {r = (n,r 2 ) : n = a + 0x, r 2 = 7 + Py, IR}, VT G T fc , 


and set 

V k := {r G (L 2 (D)) 2 : div r G L 2 (Q), rn = 0 on Tj, t\t € RT(T ) VT G 7^}, 
W k := {tp G L 2 (fl) : <p\ T e P 0 (T) VT G T k ). 
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Then the mixed finite element discretization of (2) is defined as: 


Find (Jk ,Pk) £ Vk x Wk such that for all (7k,<Pfc) £ I 4 x Wk 

[ e^JkTkdx + f pjkdiv7kdx= [ XTknds , 

J n Jn J r 0 

/ v^fcdiv Jkdx = / ftpkdx. 

Jn Jn 

The matrix associated with (3) is not coercive. To avoid this inconvenience we introduce a La- 
grange multiplier. We define 


( 3 ) 


Vi := {r € (L 2 (tt)f : At £ JET(T) VT € T fc }, 


and for £ 6 L 2 (Fo) 


A fc ^ :={//: /i e L 2 (£ k ), p\ e G P 0 (e) Ve € J (fi - £)ds = 0 Ve C r 0 }. 


Instead of (3) we now consider the mixed equilibrium discretization, 

Find (Jk,Pjt, Ajfc) G Vfe x Wk x A fc)X suc/i that for all {tv, (fik, Hk) £ Vfc x W* x Afc )0 


/ JkTkdx + VJ / pfcdiv 7kdx — / Afc7iknds= 0, 

7.CT- Jt T'cz'T. J QT 


TeT k 

/ ^Jfcdiv Jkdx 

ns-sr JT 


Ter k 


T£T k 


rcT. 


wnds 


= / /¥>*«**, 

Jn 

= 0 . 


( 4 ) 


re r fc 


As shown in [3], problem (4) has a unique solution and Jk = Jk, Pfc = Pfc holds. Moreover, 
A* is a good approximation of the solution of (2) at the interelement boundaries [2] . It is pos- 
sible to eliminate the unknowns, corresponding to Jk and Pk in the resulting system, by static 
condensation [3]. This yields a matrix (acting only on the interelement multiplier Ajt), which is a 
symmetric positive definite matrix and which is an M-matrix if the triangulation is of the weakly 
acute type (i.e. no angle > f ). 
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2.2. The Nonconforming Finite Element Formulation 


To introduce the nonconforming finite element formulation we need the following definitions: 
Let II® be the L 2 —projection from L 2 (Sk ) onto 

Afc : = {^k £ L 2 (£ k ) : p k \e £ Po(e) Ve G £ k } 

and P k be the L 2 —projection from L 2 (fl) onto 

S k := H G L 2 (fl) : v k \ T G P 0 (T) VT G T k }, 

i.e. U° k (0\e = j^j Ve G ^ and P fc( u )l T = ]^j f T udx ' VT G 

The Crouzeix-Raviart finite element space [4] is defined by 


S k := (v* 6 L 2 ((2) :v k \r € P\(T) VT eT k , v k is continuous at midpoints of edges}. 


For £ G L 2 (r 0 ) we define 


S k ,Z : = {*>* € S k : v k {m e ) = n2(£)| e> e C r 0 }. 


Notice that the standard basis functions of S k are equal to one at the midpoint of exactly 
one edge and vanish at the midpoints of all other edges. Using the arguments concerning static 
condensation in [5] , it is straightforward to prove the following lemma. 


Lemma 2.1. 


The solution X k of (4) can be written as X k = U k (w k ), where w k is the solution of the follow- 
ing nonconforming weak problem 

Find w k G Sfc iX such that for all v k G S k , o 

5^ (^fe(^)) -1 [ grad w k grad v k dx = *£(/) / (| “ 2 P°(e^ VkdX ' ^ 

m rr J T TcT. V / 


Remark 2.2. 


For w k as in Lemma 2.1. and the solution p of (2) the following error estimate [2] holds: 

_||p-Wfc||o < 7|fcfc| a (||p|ls + Plla) 

with 7 = 7(e^) independent of p and h k . 
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The Lagrange multiplier \ k is an approximation of p = e* u. In semiconductor simulations 
the range of ip is very large, so that A* is not suited for actual computations. Moreover we are 
interested in approximating the solution u of (1). Hence we introduce the following change of 
variable 

Vk = (n°(e^)) -1 A* G A fc (n o( e v.))-i x . (6) 

Denote the standard basis of S k by y? e , e E £ k . We define the linear operator E k : S k —* S k by 


EkM = n°(e^)| e ^ e Ve 6 e h . 
For / e L 2 (fl), G k (f) G L 2 (Q) is defined by 


-*(/)(§- 55^). 


Finally we arrive at the following statement: 

Lemma 2.3. 

Let C = (n° fc (e^)) 1 x G L 2 (ro). Then \i k of (6) can be written as //* = n^(u^), where u k is 
the solution of the nonconforming weak problem: 

Find u k G Skx such that for all v k G S k ,o 

Y, (- Pfc ^)) -1 f grad E k (u k ) grad v k dx = 

TeT fc Jt 

o 



Remark 2.4. 

Note that problem (7) is the usual nonconforming Crouzeix-Raviart discretization of the 
Laplace equation, if ip and / are constant on ft. 


We can use the error estimate of Remark 2.2. to obtain an estimate for the approximation 
error between the solution u k from (7) and the solution u of (1), though the result is rather un- 
satisfying. To arrive at an improved error bound, one could use the fact that two Babuska- 
Brezzi conditions hold [6] for the corresponding bilinear form. The stability and the unique solv- 
ability of the discrete problem (7) also follow. In the following we construct a multigrid algo- 
rithm for problem (7). Therefore we define the bilinear form a k on S k by 


ak(u k ,v k ) := ^2 ( p k( e ^)) 1 f grad E k (u k ) grad v k dx. 
tct, Jt 
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3. MULTIGRID METHOD 


3.1. Adaptive Mesh-Refinement Techniques 


In order to formulate the multigrid algorithm, we need a regular sequence of triangulations 
{Tk}k> o- In our refinement process, two objectives are pursued. First, in order to improve ap- 
proximation, we should refine the grid locally, where the solution behaves very badly. Second, we 
have to construct weakly acute triangulations to guarantee that the corresponding discretization 
matrix is an M-matrix. Therefore we define the strategy and rules below. Given a triangulation 
we refine its triangles as follows: 

(1) The refinement process is started by a suitable error estimator, e.g. based on residuals, 
which marks some of the triangles as red. 


(2) If a triangle is marked 

(i) red, it will be cut into four new ones by joining the midpoints of its edges, 

(ii) green, it will be cut into two new triangles by joining the midpoint of the longest edge 
to the vertex opposite to this edge, and 


(in) 


blue, it will be cut into three new triangles by joining the midpoint of its longest edge to 
the vertex opposite to this edge and to the midpoint of one of the remaining edges (see 


Fig.l) 



Figure 1. Red, 



(3) Hanging nodes are avoided using the following rules: 


(i) a triangle with three hanging nodes is marked red 

(ii) a triangle with two hanging nodes is marked blue, if one of the nodes lies on the longest 
edge of the triangle; otherwise it is marked red 

(iii) a triangle with one hanging node is marked green , if the node lies on the longest edge of 
the triangle; otherwise it is marked blue 


Note that rules (ii) and (iii) may introduce new hanging nodes. However, one can prove that 
the refinement process obeying the above rules is finite. Moreover, assuming that % has only 
isosceles right-angled triangles, then it is guaranteed that all triangulations % are weakly acute. 


6 


3.2. The Prolongation 


In order to solve problem (1), we have to find the solution Uk of the discrete problem (7). 
Since the Crouzeix-Raviart element is nonconforming and Sk-i $7 Sk-, we must construct a suit- 
able transfer operator between Sk-i and Sk- In addition, the discretization shows upwind effects 
due to the existence of strong convection in part of the domain. This also must be taken into ac- 
count. 

In [7, 8] a hierarchical basis multigrid method was used to solve a linear system arising from 
the convection diffusion equation by an upwind discretization. It was shown that the convergence 
of the hierarchical basis multigrid method depends on the strength of the convection term. When 
solving the discrete problem (7) with the multigrid algorithm [9], a similar effect can be seen in 
the numerical experiments. On the other hand, considering the one dimensional problem, one 
sees that a good interpolation has to regard the upwind effect. Therefore we introduce the fol- 
lowing weighted L 2 — projection. Define 

(«.«)* •= XI E /w<* v.cSb.ES^v (8) 

Ter fc Ter k+ 1 r 

T CT 

For all tt G PziT), T G 7^, the quadrature rule 

I udx = ® u ( m e) 

JT eC dT 

is exact, so that (8) can be written as 

(u,v) k = ^2 ( p Jfc( e ^)lf) _1 X! E k (u)(m e ) v(m e ) (9) 

TeTk TeT k+1 eC dT 

TCT 

for all ue Sk and v 6 Sk U Sk+i- 

Remark 3.1. 

Note that if v € Sk holds, (9) reduces to the equation 

(«.»)» = E( P °( e ' > )W I T E n l(e*)U u(m e ) v(m'). 

TeTfc eC dT 

Moreover, if rp is constant, we have 

(u,v) k = (u,v) V u e S k , v 6 S k USjfc + i, 
where (u,v) := J n uv dx denotes the usual L 2 — inner product. 


o 
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From (9) it follows that the standard basis functions of S k are mutually orthogonal with re- 
spect to the inner product •) k- Therefore we can obtain an easily computable prolongation op- 

erator P^-i '■ Sk - 1 ~ ► Sk by 


( Pk-\Uk-i,vk)k = (uk-i,vk)k-i V uk-i e Sk- 1 , vk e Sk- 


it is straightforward to prove the following lemma: 


Lemma 3.2. 

Let € Sfc-i, then: 

* , , w , , ng(e^)[e x -l (E k-iUk- i)(m e )lf 

If e G S k then ( P k-l u k-l)( m e) - ( pO^^| T ) Pk-il^)\f ’ 

where T (resp. f) is the triangle in % (resp. T k -i) with e C dT (resp. e C df). 
If e e £k then 

(pk v™ i - ni*i n ° (e ^ - + it r i 

(Ffc_iWfc-i)("ie) - (|T | p0^y L + I '/f(eV-)| 

ftmLi i. p 'k—l‘^'k—l)(j n e)\f L , \ 

(|T 1 

where T L ,T R (resp. T L ,T R ) are the two triangles in T k (resp. T k -\) with T L fl T R 
= e and T L C T L , T R C T R . (see Fig.2) 


Remark 3.3. 


If ^ of Lemma 3.2 is constant, we have the usual L 2 -projection as given in [9]. The 

coefficients , A: > 0, are also computed during the construction of the stiffness matrix, 

hence the interpolation is not very expensive. On the other hand, as shown in [5] the coefficients 

k> 0, introduce an upwind effect; i.e. the coefficient corresponding to the downwind 
- 

node is equal to zero. 
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Figure 2. Interpolation. 


3.3. The Smoother 


A suitable smoother for the system (7) is given in [10] by a Gauss-Seidel-iteration with de- 
coupling. This smoother is confined to special triangulations and does not allow adaptive grid 
refinements. Another candidate for problems with strong convection terms is the ILU-iteration. 
Here we restrict ourselves to a variant of the ILU-iteration. The ILU-decomposition of the linear 
system A k , related to problem (7) and the standard basis of the Crouzeix-Raviart finite element 
space Sk£, can be written as 


Ak — -Lfct/fc Dki 


where L k , Uk and Dk are given by the sparsity pattern of Ak- Denote by otk = ( a e)eef fc the coeffi- 
cient vector of Uk = a e^Pe € S k x and by bk the right hand side. Then the ILU-iteration is 

given by: 


a® an arbitrary starting vector, w € (0, 1] , 

4 = 4" 1 + uiLkUk)- 1 ^ - Afeoi" 1 ), Vi = !,•■•. 


In order to get a good smoothing rate, we must optimize the factor 


Ukjcj ~ Qfc)lb 


as mentioned in [11]. Here a* k is the solution of A k a k = b k - Therefore, by computing the optimal 
damping parameter u in every step, our final smoothing algorithm is 
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Algorithm 3.4. 


a® an arbitrary starting vector, r® = 6* — A* a®, 
for i = 1, • • • , compute: 

4 " 


vl 1 = A k d* k l , 


. 1 — l 7 *,1 — 1 


u 


i - 1 _ u k 






a* = a l fc 1 ++>' *4 *> 


rl = rl 1 - u 


,i — l„,i— 1 


end. 


o 


Remark 3.5. 

Algorithm 3.4. can be interpreted as a minimal residual method with ILTJ-preconditioning. 

o 

3.4. Multigrid Algorithm 

Now we are in the position to formulate our multigrid algorithm. 

Algorithm 3.6. (One MG-iteration at level k) 

(1) Pre-smoothing: Given u\ = Yl e e£ k Q eV ? e £ For i = 1 , • ■ • , compute using Algo- 
rithm 3.4. 

(2) Coarse-Grid Correction: Denote by u’^_ i 6 Sk- i,o the solution of the coarse grid problem 


— (^fc(/)t Ik—l V k— l) a k{ u i c 1 j -ffc— l^fc— l) — 1 £ l, 0 . 


(*) 


If k = 1, set Uk - 1 = u*k-i- If k > 1, compute an approximation Uk-i to u^__ 1 by applying 
/x = 1 or fj, = 2 iterations of the algorithm at level k — 1 to problem (*) and starting value 0. 
Set 


u 


"i+i — „,"i 


+ P k-l^k- 1- 


I 


(3) Post-smoothing: Apply iterations of Algorithm 3.4. to «^ 1+1 . 


o 
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Remark 3.7. 

So fax, there exists no convergence proof for Algorithm 3.6. The standard convergence analy- 
sis, as in [9, 11], cannot be used here, because the bilinear form ajt(., .) is unsymmetric. 


o 

4. NUMERICAL RESULTS 


In this section we present three numerical examples which demonstrate the behaviour of the 
proposed multigrid method. In all experiments we measure the performance of a method by the 
arithmetic mean of the convergence rates 



r o r r o » 

r k r k 


where r l k is the defect of the i — th iteration. 

The first model problem is taken from the papers of Brezzi, Marini and Pietra [3, 5]. We 
consider the domain fl := (0, 1) x (0, 1) with Neumann boundary 


T t := {(x,y) : ((x = 1) A {y < 0.75)) V ((y = 1) A (x < 0.75))} 


and Dirichlet boundary 

^o(x,y) 


l 


with 


Tq := dfl \ Tj, right hand side / = 0 and potential i/ 7 defined as ip(x,y) := 


V'oOuy) 


0.0 if 0.0 < r < 0.8 

r — 0.8 if 0.8 < r < 0.9 with r := \/(x — l) 2 + (y — l) 2 . 

0.1 if 0.9 < r 


On r 0 we have g(x,y ) = 0ifx = 0ory = 0 and g(x,y) = 1 otherwise. We use the initial 
triangulation % as given in Fig. 3. and refine every triangulation by marking all triangles as red 
(uniform refinement). The numerical solution for / = 10 6 and a locally refined grid is shown in 
Fig.4. 




Figure 3. Initial triangulation 1. 


Figure 4. Numerical Solution. 
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We test our multigrid algorithm 3.6. with two pre- and two post-smoothing steps {v\ = u 2 - 
2) and with different values of y, (only smoothing : n = 0; V-cycle : fx = 1; W-cycle : y, = 2) 
for problems with varying k max {k max = 1, • • • , 5). The corresponding convergence rates for l = 
10 and l = 10® are given in Tab.l and Tab.2 respectively. In all experiments we used the same 
arbitrary starting vector. 



Table 2. Convergence rates (l = 10 6 ) 


In Tab.3 we show the results for k max — 5 and with varying y and I = 10 m . 



Table 3. Convergence rates (k m ax — 5) 


In the second experiment we take 

e -V»(*.w) 

^ ^ (cosh(100 (r — 0.65))) 2 ’ 

with ip(x,y) = 10 3 (1 + tanh(100 (r - 0.65))) and r = ^/(x - l) 2 + (y - l) 2 . Again we chose 
fl = (0,1) x (0, 1). The Dirichlet boundary To = dO, and g(x,y) = (x + y) e~^ x,y y The exact 
solution is given by 

u(x,y) = {x + y) e~M x ' v) . 


The numerical solution is shown in Fig.7. We used three different coarse grids, as given in 
Fig.3, Fig.5 and Fig.6, to show that the Algorithm 3.6. does not depend on the orientation of the 
grid. For uniform refinement and k max = 5 Tab.4 shows the results with varying y, (jj, = 0, 1, 2). 
In Tab.4 we also show the results for fc max = 6 and adaptive refinement of the grid (see Fig.8). 



mm 





Figure 5. Initial triangulation 2. 


Figure 6. Initial triangulation 3. 



Figure 7. Numerical solution. Figure 8. Adaptive refined grid (k = 4). 


grid 

init. triang. 1 

init. triang. 2 

init. triang. 3 

loc. ref. . 

fi = 0 

.900 

.903 

.891 

.905 

/i= 1 

.311 

.225 

.216 

.409 

H = 2 

.157 

.087 

.128 

.355 


Table 4. Convergence rates 


Finally we consider am experiment with a real-world problem. Fig.9 shows the schematic 
structure of the doping of a thyristor. With an existing simulation program (ABBPISCES) we 
computed the solution u of (1) and the potential ip of the coupled stationary semiconductor 
equations for a blocking-state (see Fig. 11 resp. Fig. 12) and an on-state of the thyristor (see 
Fig. 13 resp. Fig. 14). The so computed potential ip was substituted into equation (1) and the re- 
sulting system was solved with our multigrid algorithm. Fig. 10 shows the grid for an adaptive 
refinement (k = 5). Finally Tab.5 shows the convergence rates for Algorithm 3.6. with a suitable 
number of pre- and post-smoothing steps, with varying fj, (/i = 0, 1,2) and k mgLX = 7. 



Figure 9. Schematic structure of the doping. 











Figure 11. Solution (log). 


Figure 12. Potential. 



Figure 13. Solution (log). Figure 14. Potential. 


state 

blocking 

on 

p = 0 

.843 

.828 

p = 1 (i/j = U2 = 22) 

.249 

.112 

p = 2 (i/i = i/ 2 = 9) 

.108 

.121 


Table 5. Convergence rates 
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